/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file
    Implementation of the density_generator classes. */

// $Id$

#include "density_generator.h"
#include "preferences.h"
#include "multiphase.h"
#include "mcrx-debug.h"

using namespace mcrx;

uniform_density_generator::uniform_density_generator (T_float dtg, 
						      T_float mtg) :
  dtg_(dtg), mtg_(mtg), use_mp_(false)
{
  assert(dtg>=0);
  assert(mtg>=0);
}


uniform_density_generator::uniform_density_generator (Preferences& p) :
  dtg_ (p.defined("dust_to_gas_ratio")?
	p.getValue("dust_to_gas_ratio", double  ()): 0),
  mtg_ (p.defined("dust_to_metal_ratio")?
	p.getValue("dust_to_metal_ratio", double  ()): 0),
  use_mp_ (p.defined("use_multiphase_model")?
	   p.getValue("use_multiphase_model", bool  ()): false),
  mp_(use_mp_?p.getValue("multiphase_t0_star", double ()):0,
      use_mp_?p.getValue("multiphase_rho_th", double ()):0)
{
  assert(dtg_>=0);
  assert(mtg_>=0);

  if(use_mp_) {
    std::cout << "Using multiphase model to generate dust densities." 
	      << std::endl;
  }
}


/** Since this type of density generator only has one type of dust,
    this function is pretty trivial.  */
T_densities uniform_density_generator::make_density_vector (T_float gd,
							   T_float md) const
{
  assert (gd >= 0);
  assert (gd == gd);
  assert (md >= 0);
  assert (md == md);

  T_densities rho (1);
  rho = dtg_*gd + mtg_*md;

  if(use_mp_) {
    // if we are using the multiphase model, we only include the
    // density of the hot phase.  This depends on the gas density.
    // 1Msun/kpc3 = 6.77e-32 g/cm3
    const double cf = mp_.cold_fraction(gd*6.77e-32).first;

    DEBUG(3, std::cout << "Gas density " << gd << " cold fraction " 
	  << cf << std::endl;);
    rho *= (1-cf);
  }

  ASSERT_ALL(rho>=0);
  return rho;
}


  
